NASA-CR-205326 


Development of a Composite Tailoring Technique for Airplane Wing 


Progress Report/Renewal Proposal 
Grant No. NAG2-908 
Technical Monitor: Dr. Hiro Miura 



NASA Ames Research Center 



Principal Investigator: Aditi Chattopadhyay 
Graduate Research Associate: Ratneshwar Jha 
Department of Mechanical and Aerospace Engineering 
Arizona State University 
Tempe, Arizina 85287-6106 



Objectives 

Composite wings have become increasingly popular in recent years due to significant 
potential of weight savings and stiffness tailoring. For aircraft wings made out of composite 
materials, aeroelastic tailoring presents an opportunity to enhance the aircraft performance by 
utilizing their unique stiffness and strength properties. Therefore, the goal of the current researh 
are as follows. 

(1) Development of a new composite beam modeling technique to represent the principal load- 
carrying member in the wing. The theory includes the effect of through-the-thickness shear 
deformation which is important in laminated composites and is not included in the classical theory. 
A refined higher-order theory is used to describe the displacement field in each wall and 
appropriate boundary conditions are imposed to ensure the satisfaction of stress-free boundary 
conditions at the free surfaces. The theory is implemented using the finite element method. 

(2) Development of a formal design optimization procedure to investigate the effect of composite 
tailoring on aeroelastic stability and structural characteristics of airplane wings. The composite 
structural analysis is coupled with unsteady aerodynamic analysis to solve the coupled aeroelastic 
equations of motion. The unsteady aerodynamic computations are performed using a panel code 
based on the Doublet Lattice Method and flutter/divergence speed is obtained using the V-g 
method. A hybrid optimization technique is implemented for the optimization to simultaneously 
include continuous and discrete design variables. 

(3) Use the developd procedure to perform design optimization studies on realistic airplane 
configurations to investigate the various aeroelastic/structural/dynamic design issues. 

Approach 

Structural Analysis 

The accurate and efficient prediction of structural response is very important for the 
investigation of aeroelastic tailoring using composite structures. The analysis of aircraft wings can 
be accomplished either through a detailed investigation of structural sections comprising spars, 
webs, ribs etc., or through the use of reduced structural models. The detailed analysis is 
computationally very expensive and is often impractical in design optimization and/or trade-off 
studies. Reduced structural models are more frequently used which include equivalent plate 
models and box beam models. Between these two, box beam models more closely represent real 
wing structures and more accurately account for elastic couplings. 
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A rectangular composite box beam model with taper and sweep is developed to represent 
the load carrying member of an aircraft wing (Fig.l). The single-celled composite box beam 
model is based on a higher-order composite laminate theory [10] and accounts for the distributions 
of shear strains through the thickness of each wall. The displacement f.eld for each wall section is 
described by bending, warping and in-plane stretching. 

For each of the individual plates, the higher-order displacement field is defined in local 
coordinate system as follows (Fig. 2 ). 

u(x,y,z,t) = UQ(x,y,t ) + zy x (x,y,t) 

+z 2 ^ x (x,y,t) + z 2 C x (x,y,t) 

v(x.y,z,t)=v 0 (x,y,t) + zy y (x,y,t) 

+z 2 $ y (x,y,t) + z 3 ^ y (x,y,t) 


where u o, v o and w o denote the displacements of a point (x, y) on the midplane and V*_and Vy 
are the rotations of the normal to the midplane about the y and x axes, respectively. The higher- 
order terms x , <=x t ^ and ; y represent beam warping in each plane. Making the assumption of 

small displacements and rotations, a linear strain-displacement relationship is used. The following 
constitutive relation is used for plates made of orthotropic materials. 
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The higher order terms are determined using the condition that the transverse shear 
s resses, a xz and a yz , vanish on the plate top and bottom surfaces. For composite plates made up 

of layers of orthotropic lamina, these cond.tions are equivalent to the requirement that the 
corresponding strains be zero on the surfaces. By making substitution for Vx _and Vy in terms of 

6 X and 0 y which are the shear angles at midplane about x and y axes respectively, [he following 
refined higher-order displacement field is obtained. 
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W = w 0 


Based on the above displacement field, a four-node plate element is developed. 

In the finite element formulation, the plate displacements at the mid-plane, u 0 ,v 0 and w 0 

are interpolated by Hemute cubic functions and the shear angles, <J» 1 and <}> y are interpolated using 
bilinear functions. For a higher-order plate element, there are 1 1 degrees of freedom at each node. 

The governing equations of motion for an individual plate is derived using Hamilton’s 
principle. 

J‘ 2 5[U -V + W nc ]dt = 0 

( 4 ) 


where U, V and W denote the kinetic energy, the strain energy and the work done by external 
forces, respectively. Using the constitutive relations along with the strain-displacement relations, 
the element stiffness matrix K e , the mass matrix M e and the forcing vector F e are derived from 
Eqn. 4 as follows. 

K e = J Ve B e T C m B e dV e 
M e =/ Ve pN e T N £ dV e 
F e = Ja N e T p( X ’ y> t)4A e 

( 5 ) 

where v e and A e represent element volume and surface area, respectively and p denotes material 
density. The matrix C m is material stiffness matrix and p is the air pressure. Matrices B e and N e 
relate the generalized coordinates to strains and displacements, respectively. 

The construction of the box beam from plate elements is shown in Fig. 3. The quantities 
u, v, and w are displacements along x, y and z axis, respectively and e x , e y and 0 2 are 

rotations along these directions. To make stiffness transformation possible, continuity of 
displacements and rotations are imposed at each of the four corners while the generalized forces 
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corresponding to higher order warping terms are set to zero. Through the use of coordinate 


transformation, the reduced stiffness matrix is expressed in the global form. Assembly of the 


element matrices leads to the following equation of motion for a general 
system. 

M x + Cx + Kx = Q\ +Q 2 


structural-aerodynamic 

(6) 


where M, C and K denote global mass, damping and stiffness matrices respectively. The vector x 
represents structural elastic deformation. The quantities Q, and Q 2 denote aeroelastic forces and 
other forces due to gust, control surface motion etc. respectively 

Aeroelastic Analysis 

For aeroelastic stability analysis, the damping C and the non-aeroelastic forces Q 2 are 
ignored. Assuming simple harmonic motion, that is, x = xe iox and Q x = Q ie ia * yields 

(—co 2 M + K)x = Q] 

(7) 


Q, can be expressed as a linear combination of x as follows. 

Q\ = q„F(ico)x 

( 8 ) 


where F(ico) is the aerodynamic influence coefficient. Substituting for Q, in Eqn (7) gives 
(-C0 2 M + K-q x F(ico))x = 0 

(9) 


Equation 9 represents an eigen value problem and the solution of 
|-<u 2 jW + K -q^Fiico ) |= 0 


( 10 ) 


provides the roots which determine the stability of the system. To solve the above problem, 
artificial damping is introduced and Eqn. 10 is rewritten as 


-co~M + (1 + ig)K - F(i<y) 


= 0 


(ID 


The solution of Eqn. 1 1 yields the variations of g and co with respect to q_ . At the flutter point, 
the damping g=0. 
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The V-g method of flutter prediction is the classical method which is widely used. In this 
method, the aerodynamic forces need to be calculated for real co only. However, the results are 
considered accurate only at the flutter point. 

Hybrid Optimization 

The inclusion of both continuous and discrete design variables significantly complicates the 
optimization problem. This is because the discrete design variables are not compatible with 
traditional gradient based optimization methods. Similarly, the continuous variable is not 
compatible with combinatorial optimization methods, such as branch and bound techniques, which 
require discrete values to operate. Therefore, a hybrid optimization technique developed by 

Chattopadhyay and Seeley which combines both types of design variables is used and is described 
next. 


The general continuous/discrete optimization problem can be mathematically stated as 
follows. 

Minimize f(o c ,<t> d ) 

Subject to <0 where j = 1 , 2 NCON 

Side constraints <j> < o < o 

c I c c u 

O d€ [$ v O v <t> % <D d J 

where f is the objective function, gj are the constraints, O c are the continuous design variables and 
<t> d are the discrete design variables which can be selected from among a set of q preselected 
values. The hybrid optimization procedure is based on Simulated Annealing (SA) where the 
design space is sampled by repeatedly perturbing the discrete design variables. At each iteration of 
the SA procedure, the objective function is minimized with respect to the continuous variables 
using a BFGS search algorithm. This significantly improves the efficiency of the hybrid algorithm 
by directing the search using the gradient information when available. The constrained problem is 
formulated using a penalty function approach. 


Results 

The development of the higher order box-beam model has been completed. Comparisons 
of the results have been made with experimental data (Chopra et al.), results from a quasi-analytical 
method (Chopra et al.), and the variational asymptotical approach (VABS, Hodges et al.). The 
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dimensions of the test beams are defined in Fig. 4. Beams with various ply lay-ups have been 
eva uated, of which only symmetric lay-up results are presented here (Figs. 5 and 6). The p, 
angles for the top and bottom walls are 45' and those for the side walls are 457-45' This land of 
lay-up exhibits bending-, orsional coupling, hence a bending load also generates beam twis, For 
the bending induced twist, results from the present model are closes, to the test data Good 

alS0 0bserved With CX P e ™'"' s f °' ■»« bending slope. Further results can be 


A simple optimization problem has been formulated to study the effects of aeroelastic 
adoring. A standard elliptic static aerodynamic load distribution has been assumed The objective 
, to nummtze the weight of the box beam which represents the structural member of an 

flutt /d- ° S ramtS are P Ced ° n the flUtter Speed and the maximum allowable stresses The 
“3": SP : d (Vf> ‘ S C ° nS,rained “ b ' - 450 knots e q uiva,en, air sp^ 

(ffiAS) a a flight condition of Mach 0.7 a, sea level. The Tsai-Wu failure criterion is imposed on 
the critical ply stresses a, the root section where material failute is most likely to occur. 

and in "* hybr ‘ d °“° n pr °“ dure “ '« Tables 1 and 2 

F,gs ' 7 - 10 ' Optimizatton results are compared with a reference design, which is selected 

M ed on engineering judgment. „ should be noted that the optimum design is independent of *e 

dial destgn due to probabilistic nature of the hybrid optimization ptocedure. The penalty function 

r:~ n r r r itera,ron ° f - -***-. zz*::z 

t^veral BFGS evaluations. Both the trial designs and the best * so far are presented Initially 
the flutter constrain, is violated which resuits in veiy tage values of the penalty function which are 

and^eTf T ^ ^ ^ 0pti ™ m Sta ‘ e * reached in less tha " 100 iterations 

found* 6 ' Pr ° CedUre “ ““ ate 250 better design could be 

Fi e g , ‘f' Si f nifiCan ‘ rCdUCU ° n in lhe wi * ht ° f 'be structural member of the wing (32% 

Fig.8, long w„h a large improvement in the flutter speed (75%. Fig.9, due to the optimization 

Since h Cmer, ° n iS Sa ' iSfied by the referenCC “ Wc " “ thc 0P' imal design (Fi. 10) 

Since the wing root chord for the reference and the optimal wings are nearly same (Table 7)' 

weigh, reduction is due to the smaller number of plies for the optima, wing. ThrLvh op ^tion 

of the stacking sequence, even a lower wall thickness provides higher flutter speed.' 

d ‘ tbe fret l ue ncies and modes for flutter show important trends. For the reference 

^ 15 3 C °“ pled b '" dins and las mode Wl,h a naIural frequency of 34 
a, .9 Hz. The fust torsion mode is the sixth mode with a natural frequency of 164 Hz. 


7 


At the flutter condition, the frequencies of the second and the sixth modes almost coalesce and the 
sixth mode also flutters at a slightly higher speed. For the optimal wing, flutter occurs for the fifth 
mode (at 74 Hz), which is the first torsion mode with the natural frequency of 160 Hz. The first 
four modes are bending or coupled bending/lag modes. Thus the optimization essentially stiffens 
the bending modes to increase the flutter speed. 

Since the aspect ratio and taper ratio are fixed for this study, smaller root chord also means 
smaller span and lower wing area. The optimal value ,of the wing root chord is very near to the 
minimum value specified to have high stiffness. This trend is expected in the absence of other 

design considerations such as wing loading (which affects landing/take-off speed, maneuverability 
etc.) and internal fuel volume. 


Future Work 

The V-g method of flutter prediction is only valid at the flutter point. Therefore, the s- 
domain method of flutter analysis, which provides both pre- and post-flutter history, will be 
adopted and will be integrated with the composite structural analysis procedure. Further details of 
the s-domain method is given in Annexure 1 . The developed analysis procedure will then be used 
within the design optimization loop to provide aeroelastically tailored wing designs. 
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Table 1 Wing Parameters. 


Reference Optimum 


Number of plies 

28 

18 

Root chord (in.) 

15.0 

15.4 

Wall thickness (in) 

0.14 

0.09 

Stacking Sequence 
top and bottom walls 

[0°/90°] 14 

[07-45 °]9 

side walls 

[457-45°] ]4 

[070°]q 


Table 2. Frequencies and Modes 


Mode Number 

Reference 

Optimum 


Mode 

Freq. 

Mode 

Freq . 



(Hz) 


(Hz) 

1 

B 

9.4 

B 

8.75 

2 

L, B 

34.1 

L, B 

34.7 

3 

B, L 

50.7 

B, L 

46.8 

4 

B, L 

1 16.9 

B 

107.7 

5 

L 

154.7 

T 

159.5 

6 

T 

163.6 

L, B 

170.1 


Flutter Point 2nd mode, 29 Hz 5th mode, 74 Hz 


Legend: B - Bending L - Lag 


T - Torsion 


X 


i 


load carrying box beam 
, A 



(a) 



(b) 

Figure 1 . Wing Geometry (a) top view (b) side view 
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Figure 3. Beam construction 
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Bending slope, 3w/3x (rad) 



d (in) 0.537 1.025 

c (in) 0.953 2.060 

Ply thickness (in) 0.005 0.005 

Wall thickness (in) 0.03 0.03 


Figure 4. Test beam dimensions 



Fig. 5 Bending slope of [45°]6 thin-walled beam under 1 Ib.bending load at tip 
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Experimental 


Fig. 6 Bendin, 


Quasi-Anal ytical 



-induced twist of [45°]6 thin-wailed beam underl lb. bending load at tip 
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Fig 7 Penalty function iteration history 
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Fig. 8 Box beam weight. 
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Fig. 9 Flutter Speed. 
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Annexure 1 


Flutter Analysis using the s-domain method 
Equation of motion for a general structural-aerodynamic system is given as 

MX + CX + KX = Q1 + Q2 (1) 

where M, C and K represent mass, damping and stiffness matrices respectively. Q1 and Q2 denote 
aeroelastic forces and other forces due to gust, control surface motion etc. For aeroelastic stability 
analysis, damping C and non-aeroelastic forces Q2 are ignored. Laplace transform of the resulting 
equation yields 

(s 2 M + K)X = Ql(s) (2) 

Ql(s) can be expressed as a linear combination of X as 

Ql(s) = q TO F(s)X (3 ) 

where F(s) can be regarded as the aerodynamic transfer function. Substituting for Ql(s) in 
Equation (2) gives 

(s 2 M + K- qoo F(s))X = 0 (4) 

This is an eigen value problem and the solution of 

s 2 M + K - q M F(s) = 0 (5) 

gives the roots which determine the stability of the system. For stability, all real roots should be 
negative. At flutter condition, one of the roots is purely imaginary. The main difficulty in solving 
equation (5) arises in obtaining the aerodynamic transfer function F(s). It is assumed that F(s) 
equals F(iw), where w is a real value (as for the V-g method). F(iw) is obtained from the Doublet 
Lattice Method code mentioned earlier. 

F(s) is expressed in Pade Approximant form 

F(s) = a + bs + cs 2 + Z p— (6) 

and the coefficients are evaluated using a least square fit. Then equation (5) is solved in complex s- 
domain to obtain the roots of the system. As mentioned earlier, the solution (frequency and 
damping) is valid at all speeds, not just the flutter speed, unlike the V-g method. This method is 
also advantageous in active control of aeroelastic systems. 


1 5 



